# Function to generate the record count matrix for a single hospital
generate_record_count <- function(data) {
counts <- table(data[, "n_hi"])
result_matrix <- cbind(as.numeric(names(counts)), as.numeric(counts))
colnames(result_matrix) <- c("n_hi", "frequency")
return(result_matrix)
}
# Function to generate Z_hv matrix for a single hospital
generate_Zhv_matrix <- function(data) {
record_count_matrix <- generate_record_count(data)
diagonal_blocks <- list()
for (i in 1:nrow(record_count_matrix)) {
n_hi <- record_count_matrix[i, "n_hi"]
frequency <- record_count_matrix[i, "frequency"]
identity_block <- diag(frequency)
ones_vector <- matrix(1, nrow = n_hi, ncol = 1)
kronecker_product <- kronecker(identity_block, ones_vector)
diagonal_blocks[[i]] <- kronecker_product
}
big_matrix <- do.call(Matrix::bdiag, diagonal_blocks)
big_matrix <- cbind(1, big_matrix)
return(as.matrix(big_matrix))
}
# Function to get summary stats from each site for distributed LMM
lmm.get.summary3 <- function(Y = NULL, X = NULL, Z = NULL, id.site = NULL, weights = NULL) {
if (is.null(weights)) weights <- rep(1, length(Y))
X <- as.matrix(X)
id.site <- as.character(id.site)
id.site.uniq <- unique(id.site)
px <- ncol(X)
ShXYZ <- list()
for (h in seq_along(id.site.uniq)) {
sh <- id.site.uniq[h]
wth <- weights[id.site == sh]
Xh <- X[id.site == sh, ]
Yh <- Y[id.site == sh]
Zh <- Z[h][[1]]
ShX <- t(Xh * wth) %*% Xh
ShXZ <- t(Xh * wth) %*% Zh
ShXY <- t(Xh * wth) %*% Yh
ShZ <- t(Zh * wth) %*% Zh
ShZY <- t(Zh * wth) %*% Yh
ShY <- sum(Yh^2 * wth)
Nh <- sum(id.site == sh)
ShXYZ[[sh]] <- list(ShX = ShX, ShXZ = ShXZ, ShXY = ShXY,
ShZ = ShZ, ShZY = ShZY, ShY = ShY, Nh = Nh)
}
return(ShXYZ)
}
##############################
# different site different patients
# Parameters
H <- 5 # of sites
m_hosp <- sample(100:150, H, replace = T) # of patients (1k to 3k later)
px <- 9 # of covariates
p_bin <- 5 # of binary X
p_cont <- px - p_bin # of continuous X
beta0 = 5 # intercept
beta <- c(2, 4, 6, 8, 10, 3, 5, 7, 9) # Fixed effects for covariates
sigma_e <- 1 # error variance
sigma_u <- 3 # site-level variance
sigma_v_hosp <- runif(H, min = 1, max = 5) # Varying sigma_v by hospital
# Generate data
nn <- rep(m_hosp, times = 1) # Number of patients per hospital
id.hosp <- rep(1:H, times = m_hosp) # Hospital ID
id.pat <- sequence(nn) # Patient ID
n_visits <- sample(1:30, sum(nn), replace = TRUE) # number of visits of patients
# Expand hospital and patient IDs for visits
id.visit <- sequence(n_visits)
id.hosp.expanded <- rep(id.hosp, times = n_visits)
id.pat.expanded <- rep(id.pat, times = n_visits)
# Random effects
u_h <- rnorm(H, mean = 0, sd = sigma_u) # Hospital effects
v_hi <- rnorm(sum(nn), mean = 0, sd = rep(sigma_v_hosp, times = m_hosp)) # Patient effects varying by hospital
# Expansion of hospital effects
u_h_patient <- rep(u_h, times = m_hosp)
u_h_expanded <- rep(u_h_patient, times = n_visits)
# Expansion of patient effects
v_hi_expanded <- rep(v_hi, times = n_visits)
# covariates
X_bin <- matrix(rbinom(sum(n_visits) * p_bin, size = 1, prob = 0.3), nrow = sum(n_visits), ncol = p_bin)
X_cont <- matrix(rnorm(sum(n_visits) * p_cont, mean = 0, sd = 1), nrow = sum(n_visits), ncol = p_cont)
X_hij <- cbind(X_bin, X_cont) # Combine binary & continuous covariates
epsilon_hij <- rnorm(sum(n_visits), mean = 0, sd = sigma_e)
# Compute outcome
y_hij <- beta0 + X_hij %*% beta + u_h_expanded + v_hi_expanded + epsilon_hij
three_lvl_dat <- data.table(
site = id.hosp.expanded,
patient = id.pat.expanded,
visit = id.visit,
X_hij,
Y = y_hij
) %>% data.frame()
setnames(three_lvl_dat, c("site", "patient", "visit", paste0("X", 1:px), "Y"))
# Preprocessing
# Calculate the number of visits per patient per site
visit_count <- three_lvl_dat %>%
dplyr::group_by(site, patient) %>%
dplyr::summarise(total_visits = n(), .groups = "drop")
# Reorder data (as a preprocess probably)
rearranged_data <- merge(three_lvl_dat, visit_count, by = c("site", "patient")) %>%
arrange(site, total_visits, patient) %>%
mutate(site = factor(site))
# XYZ
Y <- rearranged_data$Y
X <- as.matrix(rearranged_data[, paste0("X", 1:px)])
X <- cbind(1, X)
Z <- list()
for(i in 1:H){
count_mat = rearranged_data %>%
filter(site == i) %>%
group_by(site, patient) %>%
dplyr::summarise(n_hi = n(), .groups = 'drop')
Z[[i]] <- (generate_Zhv_matrix(count_mat))
}
id.site <- rearranged_data$site
ShXYZ <- lmm.get.summary3(Y, X, Z, id.site)
source("DLMM_engine3RI.R")
bench::mark(lmm.fit3(Y = NULL, X = NULL, Z = NULL,
id.site = NULL, weights = NULL,
pooled = F, reml = T,
common.s2 = T,
ShXYZ = ShXYZ,
corstr = 'independence',
mypar.init = NULL))
default mypar.init (var comp) = 1 1 1 1 1 1
Normal exit from bobyqa The number of function evaluations used is 1739default mypar.init (var comp) = 1 1 1 1 1 1
Normal exit from bobyqa The number of function evaluations used is 1739
Warning: Some expressions had a GC in every iteration; so filtering is disabled.
fit03.dlmm = lmm.fit3(Y = NULL, X = NULL, Z = NULL,
id.site = NULL, weights = NULL,
pooled = F, reml = T,
common.s2 = T,
ShXYZ = ShXYZ, # only need summary stats
corstr = 'independence',
mypar.init = NULL)
default mypar.init (var comp) = 1 1 1 1 1 1
Normal exit from bobyqa The number of function evaluations used is 1739
fit03.dlmm$b
[,1]
2.427526
X1 2.030760
X2 4.025651
X3 6.009985
X4 8.027462
X5 9.964476
X6 3.002865
X7 5.003628
X8 7.003889
X9 8.990027
sqrt(fit03.dlmm$V)
[1] 2.006603 2.837770 3.293975 2.701197 2.025338 2.754136
sqrt(fit03.dlmm$s2)
[1] 0.9975977
true_beta = c(beta0, beta)
true_beta
[1] 5 2 4 6 8 10 3 5 7 9
true_sigma = c(sigma_u, sigma_v_hosp)
true_sigma
[1] 3.000000 2.733776 3.263261 2.339320 2.124376 2.726359
c(sigma_e)
[1] 1
plot(fit03.dlmm$b, true_beta)
points(fit03.dlmm$b[1], true_beta[1], col = "red")
abline(a = 0, b = 1, col = "blue", lwd = 2)

plot(sqrt(fit03.dlmm$V), true_sigma)
points(sqrt(fit03.dlmm$V)[1], true_sigma[1], col = "red")
abline(a = 0, b = 1, col = "blue", lwd = 2)

source("DLMM_Engine.R")
bench::mark(lmm.fit3(Y = NULL, X = NULL, Z = NULL,
id.site = NULL, weights = NULL,
pooled = F, reml = T,
common.s2 = T,
ShXYZ = ShXYZ,
corstr = 'independence',
mypar.init = NULL))
Default mypar.init (var comp) = 1 1 1 1 1 1
Normal exit from bobyqa The number of function evaluations used is 1182
Default mypar.init (var comp) = 1 1 1 1 1 1
Normal exit from bobyqa The number of function evaluations used is 1182
Warning: Some expressions had a GC in every iteration; so filtering is disabled.
fit03.dlmm = lmm.fit3(Y = NULL, X = NULL, Z = NULL,
id.site = NULL, weights = NULL,
pooled = F, reml = T,
common.s2 = T,
ShXYZ = ShXYZ, # only need summary stats
corstr = 'independence',
mypar.init = NULL)
Default mypar.init (var comp) = 1 1 1 1 1 1
Normal exit from bobyqa The number of function evaluations used is 1182
fit03.dlmm$b
[,1]
2.427516
X1 2.030760
X2 4.025651
X3 6.009985
X4 8.027462
X5 9.964476
X6 3.002865
X7 5.003628
X8 7.003889
X9 8.990027
sqrt(fit03.dlmm$V)
[1] 2.003896 2.837504 3.294069 2.701259 2.025409 2.753977
sqrt(fit03.dlmm$s2)
[1] 0.9975991
true_beta = c(beta0, beta)
true_beta
[1] 5 2 4 6 8 10 3 5 7 9
true_sigma = c(sigma_u, sigma_v_hosp)
true_sigma
[1] 3.000000 2.733776 3.263261 2.339320 2.124376 2.726359
c(sigma_e)
[1] 1
plot(fit03.dlmm$b, true_beta)
points(fit03.dlmm$b[1], true_beta[1], col = "red")
abline(a = 0, b = 1, col = "blue", lwd = 2)

plot(sqrt(fit03.dlmm$V), true_sigma)
points(sqrt(fit03.dlmm$V)[1], true_sigma[1], col = "red")
abline(a = 0, b = 1, col = "blue", lwd = 2)

source("DLMM_Engine_Efficient.R")
bench::mark(lmm.fit3(Y = NULL, X = NULL, Z = NULL,
id.site = NULL, weights = NULL,
pooled = F, reml = T,
common.s2 = T,
ShXYZ = ShXYZ,
corstr = 'independence',
mypar.init = NULL))
Default mypar.init (var comp) = 1 1 1 1 1 1
A <- matrix(runif(1024^2), nrow = 1024)
B <- matrix(runif(1024^2), nrow = 1024)
size = 1024
bench::mark(logdet <- log(det(diag(1, size) + t(A) %*% A %*% B)))
bench::mark(logdet <- as.numeric(determinant(diag(1, size) + t(A) %*% A %*% B, logarithm = TRUE)$modulus))
LS0tCnRpdGxlOiAnQ29kZXMnCm91dHB1dDoKICBodG1sX25vdGVib29rOgogICAgY29kZV9mb2xkaW5nOiBub25lCiAgICBoaWdobGlnaHQ6IHRleHRtYXRlCiAgaHRtbF9kb2N1bWVudDoKICAgIGRmX3ByaW50OiBwYWdlZAotLS0KCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBULCBldmFsID0gVCwgY2FjaGUgPSBGLCB3YXJuaW5nID0gVCwgbWVzc2FnZSA9IFQpCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGxtZTQpCmxpYnJhcnkoTUFTUykgIApsaWJyYXJ5KGRhdGEudGFibGUpCmxpYnJhcnkoZ3JpZEV4dHJhKQpnZ3Bsb3QyOjp0aGVtZV9zZXQoZ2dwbG90Mjo6dGhlbWVfYncoKSkKa25pdHI6Om9wdHNfY2h1bmskc2V0KG91dC53aWR0aCA9ICIxMDAlIiwgZmlnLmFsaWduID0gJ2NlbnRlcicpCnNldC5zZWVkKDE5MTgpCmBgYAoKCmBgYHtyfQojIEZ1bmN0aW9uIHRvIGdlbmVyYXRlIHRoZSByZWNvcmQgY291bnQgbWF0cml4IGZvciBhIHNpbmdsZSBob3NwaXRhbApnZW5lcmF0ZV9yZWNvcmRfY291bnQgPC0gZnVuY3Rpb24oZGF0YSkgewogIGNvdW50cyA8LSB0YWJsZShkYXRhWywgIm5faGkiXSkKICByZXN1bHRfbWF0cml4IDwtIGNiaW5kKGFzLm51bWVyaWMobmFtZXMoY291bnRzKSksIGFzLm51bWVyaWMoY291bnRzKSkKICBjb2xuYW1lcyhyZXN1bHRfbWF0cml4KSA8LSBjKCJuX2hpIiwgImZyZXF1ZW5jeSIpCiAgcmV0dXJuKHJlc3VsdF9tYXRyaXgpCn0KCiMgRnVuY3Rpb24gdG8gZ2VuZXJhdGUgWl9odiBtYXRyaXggZm9yIGEgc2luZ2xlIGhvc3BpdGFsCmdlbmVyYXRlX1podl9tYXRyaXggPC0gZnVuY3Rpb24oZGF0YSkgewogIHJlY29yZF9jb3VudF9tYXRyaXggPC0gZ2VuZXJhdGVfcmVjb3JkX2NvdW50KGRhdGEpCiAgZGlhZ29uYWxfYmxvY2tzIDwtIGxpc3QoKQoKICBmb3IgKGkgaW4gMTpucm93KHJlY29yZF9jb3VudF9tYXRyaXgpKSB7CiAgICBuX2hpIDwtIHJlY29yZF9jb3VudF9tYXRyaXhbaSwgIm5faGkiXQogICAgZnJlcXVlbmN5IDwtIHJlY29yZF9jb3VudF9tYXRyaXhbaSwgImZyZXF1ZW5jeSJdCgogICAgaWRlbnRpdHlfYmxvY2sgPC0gZGlhZyhmcmVxdWVuY3kpCiAgICBvbmVzX3ZlY3RvciA8LSBtYXRyaXgoMSwgbnJvdyA9IG5faGksIG5jb2wgPSAxKQoKICAgIGtyb25lY2tlcl9wcm9kdWN0IDwtIGtyb25lY2tlcihpZGVudGl0eV9ibG9jaywgb25lc192ZWN0b3IpCiAgICBkaWFnb25hbF9ibG9ja3NbW2ldXSA8LSBrcm9uZWNrZXJfcHJvZHVjdAogIH0KCiAgYmlnX21hdHJpeCA8LSBkby5jYWxsKE1hdHJpeDo6YmRpYWcsIGRpYWdvbmFsX2Jsb2NrcykKICBiaWdfbWF0cml4IDwtIGNiaW5kKDEsIGJpZ19tYXRyaXgpCiAgcmV0dXJuKGFzLm1hdHJpeChiaWdfbWF0cml4KSkKfQoKIyBGdW5jdGlvbiB0byBnZXQgc3VtbWFyeSBzdGF0cyBmcm9tIGVhY2ggc2l0ZSBmb3IgZGlzdHJpYnV0ZWQgTE1NCmxtbS5nZXQuc3VtbWFyeTMgPC0gZnVuY3Rpb24oWSA9IE5VTEwsIFggPSBOVUxMLCBaID0gTlVMTCwgaWQuc2l0ZSA9IE5VTEwsIHdlaWdodHMgPSBOVUxMKSB7CiAgaWYgKGlzLm51bGwod2VpZ2h0cykpIHdlaWdodHMgPC0gcmVwKDEsIGxlbmd0aChZKSkKICBYIDwtIGFzLm1hdHJpeChYKQogIGlkLnNpdGUgPC0gYXMuY2hhcmFjdGVyKGlkLnNpdGUpCiAgaWQuc2l0ZS51bmlxIDwtIHVuaXF1ZShpZC5zaXRlKQogIHB4IDwtIG5jb2woWCkKCiAgU2hYWVogPC0gbGlzdCgpCiAgZm9yIChoIGluIHNlcV9hbG9uZyhpZC5zaXRlLnVuaXEpKSB7CiAgICBzaCA8LSBpZC5zaXRlLnVuaXFbaF0KICAgIHd0aCA8LSB3ZWlnaHRzW2lkLnNpdGUgPT0gc2hdCiAgICBYaCA8LSBYW2lkLnNpdGUgPT0gc2gsIF0KICAgIFloIDwtIFlbaWQuc2l0ZSA9PSBzaF0KICAgIFpoIDwtIFpbaF1bWzFdXQoKICAgIFNoWCAgPC0gdChYaCAqIHd0aCkgJSolIFhoCiAgICBTaFhaIDwtIHQoWGggKiB3dGgpICUqJSBaaAogICAgU2hYWSA8LSB0KFhoICogd3RoKSAlKiUgWWgKICAgIFNoWiAgPC0gdChaaCAqIHd0aCkgJSolIFpoCiAgICBTaFpZIDwtIHQoWmggKiB3dGgpICUqJSBZaAogICAgU2hZICA8LSBzdW0oWWheMiAqIHd0aCkKICAgIE5oIDwtIHN1bShpZC5zaXRlID09IHNoKQoKICAgIFNoWFlaW1tzaF1dIDwtIGxpc3QoU2hYID0gU2hYLCBTaFhaID0gU2hYWiwgU2hYWSA9IFNoWFksCiAgICAgICAgICAgICAgICAgICAgICAgIFNoWiA9IFNoWiwgU2haWSA9IFNoWlksIFNoWSA9IFNoWSwgTmggPSBOaCkKICB9CgogIHJldHVybihTaFhZWikKfQpgYGAKCgpgYGB7cn0KIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiMgZGlmZmVyZW50IHNpdGUgZGlmZmVyZW50IHBhdGllbnRzCgojIFBhcmFtZXRlcnMKSCA8LSA1ICAjIG9mIHNpdGVzCm1faG9zcCA8LSBzYW1wbGUoNTA6MTAwLCBILCByZXBsYWNlID0gVCkgIyBvZiBwYXRpZW50cyAoMWsgdG8gM2sgbGF0ZXIpCgoKcHggPC0gOSAgIyBvZiBjb3ZhcmlhdGVzCnBfYmluIDwtIDUgICMgb2YgYmluYXJ5IFgKcF9jb250IDwtIHB4IC0gcF9iaW4gICMgb2YgY29udGludW91cyBYCgoKYmV0YTAgPSA1ICMgaW50ZXJjZXB0CmJldGEgPC0gYygyLCA0LCA2LCA4LCAxMCwgMywgNSwgNywgOSkgICMgRml4ZWQgZWZmZWN0cyBmb3IgY292YXJpYXRlcwoKc2lnbWFfZSA8LSAxICAjIGVycm9yIHZhcmlhbmNlCnNpZ21hX3UgPC0gMyAjIHNpdGUtbGV2ZWwgdmFyaWFuY2UKc2lnbWFfdl9ob3NwIDwtIHJ1bmlmKEgsIG1pbiA9IDEsIG1heCA9IDUpICAjIFZhcnlpbmcgc2lnbWFfdiBieSBob3NwaXRhbAoKCiMgR2VuZXJhdGUgZGF0YQpubiA8LSByZXAobV9ob3NwLCB0aW1lcyA9IDEpICAjIE51bWJlciBvZiBwYXRpZW50cyBwZXIgaG9zcGl0YWwKaWQuaG9zcCA8LSByZXAoMTpILCB0aW1lcyA9IG1faG9zcCkgICAjIEhvc3BpdGFsIElECmlkLnBhdCA8LSBzZXF1ZW5jZShubikgICAgICAgICAgICAgICAgIyBQYXRpZW50IElECm5fdmlzaXRzIDwtIHNhbXBsZSgxOjMwLCBzdW0obm4pLCByZXBsYWNlID0gVFJVRSkgICMgbnVtYmVyIG9mIHZpc2l0cyBvZiBwYXRpZW50cwoKIyBFeHBhbmQgaG9zcGl0YWwgYW5kIHBhdGllbnQgSURzIGZvciB2aXNpdHMKaWQudmlzaXQgPC0gc2VxdWVuY2Uobl92aXNpdHMpCmlkLmhvc3AuZXhwYW5kZWQgPC0gcmVwKGlkLmhvc3AsIHRpbWVzID0gbl92aXNpdHMpCmlkLnBhdC5leHBhbmRlZCA8LSByZXAoaWQucGF0LCB0aW1lcyA9IG5fdmlzaXRzKQoKIyBSYW5kb20gZWZmZWN0cwp1X2ggPC0gcm5vcm0oSCwgbWVhbiA9IDAsIHNkID0gc2lnbWFfdSkgICMgSG9zcGl0YWwgZWZmZWN0cwp2X2hpIDwtIHJub3JtKHN1bShubiksIG1lYW4gPSAwLCBzZCA9IHJlcChzaWdtYV92X2hvc3AsIHRpbWVzID0gbV9ob3NwKSkgICMgUGF0aWVudCBlZmZlY3RzIHZhcnlpbmcgYnkgaG9zcGl0YWwKCiMgRXhwYW5zaW9uIG9mIGhvc3BpdGFsIGVmZmVjdHMKdV9oX3BhdGllbnQgPC0gcmVwKHVfaCwgdGltZXMgPSBtX2hvc3ApCnVfaF9leHBhbmRlZCA8LSByZXAodV9oX3BhdGllbnQsIHRpbWVzID0gbl92aXNpdHMpCgojIEV4cGFuc2lvbiBvZiBwYXRpZW50IGVmZmVjdHMKdl9oaV9leHBhbmRlZCA8LSByZXAodl9oaSwgdGltZXMgPSBuX3Zpc2l0cykKCiMgY292YXJpYXRlcwpYX2JpbiA8LSBtYXRyaXgocmJpbm9tKHN1bShuX3Zpc2l0cykgKiBwX2Jpbiwgc2l6ZSA9IDEsIHByb2IgPSAwLjMpLCBucm93ID0gc3VtKG5fdmlzaXRzKSwgbmNvbCA9IHBfYmluKQpYX2NvbnQgPC0gbWF0cml4KHJub3JtKHN1bShuX3Zpc2l0cykgKiBwX2NvbnQsIG1lYW4gPSAwLCBzZCA9IDEpLCBucm93ID0gc3VtKG5fdmlzaXRzKSwgbmNvbCA9IHBfY29udCkKWF9oaWogPC0gY2JpbmQoWF9iaW4sIFhfY29udCkgICMgQ29tYmluZSBiaW5hcnkgJiBjb250aW51b3VzIGNvdmFyaWF0ZXMKCmVwc2lsb25faGlqIDwtIHJub3JtKHN1bShuX3Zpc2l0cyksIG1lYW4gPSAwLCBzZCA9IHNpZ21hX2UpCgojIENvbXB1dGUgb3V0Y29tZQp5X2hpaiA8LSBiZXRhMCArIFhfaGlqICUqJSBiZXRhICsgdV9oX2V4cGFuZGVkICsgdl9oaV9leHBhbmRlZCArIGVwc2lsb25faGlqCgp0aHJlZV9sdmxfZGF0IDwtIGRhdGEudGFibGUoCiAgc2l0ZSA9IGlkLmhvc3AuZXhwYW5kZWQsCiAgcGF0aWVudCA9IGlkLnBhdC5leHBhbmRlZCwKICB2aXNpdCA9IGlkLnZpc2l0LAogIFhfaGlqLAogIFkgPSB5X2hpagopICU+JSBkYXRhLmZyYW1lKCkKCnNldG5hbWVzKHRocmVlX2x2bF9kYXQsIGMoInNpdGUiLCAicGF0aWVudCIsICJ2aXNpdCIsIHBhc3RlMCgiWCIsIDE6cHgpLCAiWSIpKQoKCiMgUHJlcHJvY2Vzc2luZwojIENhbGN1bGF0ZSB0aGUgbnVtYmVyIG9mIHZpc2l0cyBwZXIgcGF0aWVudCBwZXIgc2l0ZQp2aXNpdF9jb3VudCA8LSB0aHJlZV9sdmxfZGF0ICU+JQogIGRwbHlyOjpncm91cF9ieShzaXRlLCBwYXRpZW50KSAlPiUKICBkcGx5cjo6c3VtbWFyaXNlKHRvdGFsX3Zpc2l0cyA9IG4oKSwgLmdyb3VwcyA9ICJkcm9wIikKCiMgUmVvcmRlciBkYXRhIChhcyBhIHByZXByb2Nlc3MgcHJvYmFibHkpCnJlYXJyYW5nZWRfZGF0YSA8LSBtZXJnZSh0aHJlZV9sdmxfZGF0LCB2aXNpdF9jb3VudCwgYnkgPSBjKCJzaXRlIiwgInBhdGllbnQiKSkgJT4lCiAgYXJyYW5nZShzaXRlLCB0b3RhbF92aXNpdHMsIHBhdGllbnQpICU+JQogIG11dGF0ZShzaXRlID0gZmFjdG9yKHNpdGUpKQoKCiMgWFlaClkgPC0gcmVhcnJhbmdlZF9kYXRhJFkKWCA8LSBhcy5tYXRyaXgocmVhcnJhbmdlZF9kYXRhWywgcGFzdGUwKCJYIiwgMTpweCldKQpYIDwtIGNiaW5kKDEsIFgpClogPC0gbGlzdCgpCgpmb3IoaSBpbiAxOkgpewogIGNvdW50X21hdCA9IHJlYXJyYW5nZWRfZGF0YSAlPiUKICAgICAgICAgICAgICAgICAgICBmaWx0ZXIoc2l0ZSA9PSBpKSAlPiUKICAgICAgICAgICAgICAgICAgICBncm91cF9ieShzaXRlLCBwYXRpZW50KSAlPiUKICAgICAgICAgICAgICAgICAgICBkcGx5cjo6c3VtbWFyaXNlKG5faGkgPSBuKCksIC5ncm91cHMgPSAnZHJvcCcpCgogIFpbW2ldXSA8LSAoZ2VuZXJhdGVfWmh2X21hdHJpeChjb3VudF9tYXQpKQp9CgppZC5zaXRlIDwtIHJlYXJyYW5nZWRfZGF0YSRzaXRlCgpTaFhZWiA8LSBsbW0uZ2V0LnN1bW1hcnkzKFksIFgsIFosIGlkLnNpdGUpCmBgYAoKCmBgYHtyfQpzb3VyY2UoIkRMTU1fZW5naW5lM1JJLlIiKQpiZW5jaDo6bWFyayhsbW0uZml0MyhZID0gTlVMTCwgWCA9IE5VTEwsIFogPSBOVUxMLAogICAgICAgICAgICAgICAgICAgICBpZC5zaXRlID0gTlVMTCwgd2VpZ2h0cyA9IE5VTEwsCiAgICAgICAgICAgICAgICAgICAgIHBvb2xlZCA9IEYsIHJlbWwgPSBULAogICAgICAgICAgICAgICAgICAgICBjb21tb24uczIgPSBULAogICAgICAgICAgICAgICAgICAgICBTaFhZWiA9IFNoWFlaLAogICAgICAgICAgICAgICAgICAgICBjb3JzdHIgPSAnaW5kZXBlbmRlbmNlJywKICAgICAgICAgICAgICAgICAgICAgbXlwYXIuaW5pdCA9IE5VTEwpKQoKZml0MDMuZGxtbSA9IGxtbS5maXQzKFkgPSBOVUxMLCBYID0gTlVMTCwgWiA9IE5VTEwsCiAgICAgICAgICAgICAgICAgICAgICBpZC5zaXRlID0gTlVMTCwgd2VpZ2h0cyA9IE5VTEwsCiAgICAgICAgICAgICAgICAgICAgICBwb29sZWQgPSBGLCByZW1sID0gVCwKICAgICAgICAgICAgICAgICAgICAgIGNvbW1vbi5zMiA9IFQsCiAgICAgICAgICAgICAgICAgICAgICBTaFhZWiA9IFNoWFlaLCAgIyBvbmx5IG5lZWQgc3VtbWFyeSBzdGF0cwogICAgICAgICAgICAgICAgICAgICAgY29yc3RyID0gJ2luZGVwZW5kZW5jZScsCiAgICAgICAgICAgICAgICAgICAgICBteXBhci5pbml0ID0gTlVMTCkKCgpmaXQwMy5kbG1tJGIKc3FydChmaXQwMy5kbG1tJFYpCnNxcnQoZml0MDMuZGxtbSRzMikKCgoKdHJ1ZV9iZXRhID0gYyhiZXRhMCwgYmV0YSkKdHJ1ZV9iZXRhCnRydWVfc2lnbWEgPSBjKHNpZ21hX3UsIHNpZ21hX3ZfaG9zcCkKdHJ1ZV9zaWdtYQpjKHNpZ21hX2UpCgoKcGxvdChmaXQwMy5kbG1tJGIsIHRydWVfYmV0YSkKcG9pbnRzKGZpdDAzLmRsbW0kYlsxXSwgdHJ1ZV9iZXRhWzFdLCBjb2wgPSAicmVkIikKYWJsaW5lKGEgPSAwLCBiID0gMSwgY29sID0gImJsdWUiLCBsd2QgPSAyKQoKcGxvdChzcXJ0KGZpdDAzLmRsbW0kViksIHRydWVfc2lnbWEpCnBvaW50cyhzcXJ0KGZpdDAzLmRsbW0kVilbMV0sIHRydWVfc2lnbWFbMV0sIGNvbCA9ICJyZWQiKQphYmxpbmUoYSA9IDAsIGIgPSAxLCBjb2wgPSAiYmx1ZSIsIGx3ZCA9IDIpCmBgYAoKCmBgYHtyfQpzb3VyY2UoIkRMTU1fRW5naW5lLlIiKQpiZW5jaDo6bWFyayhsbW0uZml0MyhZID0gTlVMTCwgWCA9IE5VTEwsIFogPSBOVUxMLAogICAgICAgICAgICAgICAgICAgICBpZC5zaXRlID0gTlVMTCwgd2VpZ2h0cyA9IE5VTEwsCiAgICAgICAgICAgICAgICAgICAgIHBvb2xlZCA9IEYsIHJlbWwgPSBULAogICAgICAgICAgICAgICAgICAgICBjb21tb24uczIgPSBULAogICAgICAgICAgICAgICAgICAgICBTaFhZWiA9IFNoWFlaLAogICAgICAgICAgICAgICAgICAgICBjb3JzdHIgPSAnaW5kZXBlbmRlbmNlJywKICAgICAgICAgICAgICAgICAgICAgbXlwYXIuaW5pdCA9IE5VTEwpKQoKCmZpdDAzLmRsbW0gPSBsbW0uZml0MyhZID0gTlVMTCwgWCA9IE5VTEwsIFogPSBOVUxMLAogICAgICAgICAgICAgICAgICAgICAgaWQuc2l0ZSA9IE5VTEwsIHdlaWdodHMgPSBOVUxMLAogICAgICAgICAgICAgICAgICAgICAgcG9vbGVkID0gRiwgcmVtbCA9IFQsCiAgICAgICAgICAgICAgICAgICAgICBjb21tb24uczIgPSBULAogICAgICAgICAgICAgICAgICAgICAgU2hYWVogPSBTaFhZWiwgICMgb25seSBuZWVkIHN1bW1hcnkgc3RhdHMKICAgICAgICAgICAgICAgICAgICAgIGNvcnN0ciA9ICdpbmRlcGVuZGVuY2UnLAogICAgICAgICAgICAgICAgICAgICAgbXlwYXIuaW5pdCA9IE5VTEwpCgoKZml0MDMuZGxtbSRiCnNxcnQoZml0MDMuZGxtbSRWKQpzcXJ0KGZpdDAzLmRsbW0kczIpCgoKCnRydWVfYmV0YSA9IGMoYmV0YTAsIGJldGEpCnRydWVfYmV0YQp0cnVlX3NpZ21hID0gYyhzaWdtYV91LCBzaWdtYV92X2hvc3ApCnRydWVfc2lnbWEKYyhzaWdtYV9lKQoKCnBsb3QoZml0MDMuZGxtbSRiLCB0cnVlX2JldGEpCnBvaW50cyhmaXQwMy5kbG1tJGJbMV0sIHRydWVfYmV0YVsxXSwgY29sID0gInJlZCIpCmFibGluZShhID0gMCwgYiA9IDEsIGNvbCA9ICJibHVlIiwgbHdkID0gMikKCnBsb3Qoc3FydChmaXQwMy5kbG1tJFYpLCB0cnVlX3NpZ21hKQpwb2ludHMoc3FydChmaXQwMy5kbG1tJFYpWzFdLCB0cnVlX3NpZ21hWzFdLCBjb2wgPSAicmVkIikKYWJsaW5lKGEgPSAwLCBiID0gMSwgY29sID0gImJsdWUiLCBsd2QgPSAyKQpgYGAKCgoKCmBgYHtyfQpzb3VyY2UoIkRMTU1fRW5naW5lX0VmZmljaWVudC5SIikKYmVuY2g6Om1hcmsobG1tLmZpdDMoWSA9IE5VTEwsIFggPSBOVUxMLCBaID0gTlVMTCwKICAgICAgICAgICAgICAgICAgICAgaWQuc2l0ZSA9IE5VTEwsIHdlaWdodHMgPSBOVUxMLAogICAgICAgICAgICAgICAgICAgICBwb29sZWQgPSBGLCByZW1sID0gVCwKICAgICAgICAgICAgICAgICAgICAgY29tbW9uLnMyID0gVCwKICAgICAgICAgICAgICAgICAgICAgU2hYWVogPSBTaFhZWiwKICAgICAgICAgICAgICAgICAgICAgY29yc3RyID0gJ2luZGVwZW5kZW5jZScsCiAgICAgICAgICAgICAgICAgICAgIG15cGFyLmluaXQgPSBOVUxMKSkKCgpmaXQwMy5kbG1tID0gbG1tLmZpdDMoWSA9IE5VTEwsIFggPSBOVUxMLCBaID0gTlVMTCwKICAgICAgICAgICAgICAgICAgICAgIGlkLnNpdGUgPSBOVUxMLCB3ZWlnaHRzID0gTlVMTCwKICAgICAgICAgICAgICAgICAgICAgIHBvb2xlZCA9IEYsIHJlbWwgPSBULAogICAgICAgICAgICAgICAgICAgICAgY29tbW9uLnMyID0gVCwKICAgICAgICAgICAgICAgICAgICAgIFNoWFlaID0gU2hYWVosICAjIG9ubHkgbmVlZCBzdW1tYXJ5IHN0YXRzCiAgICAgICAgICAgICAgICAgICAgICBjb3JzdHIgPSAnaW5kZXBlbmRlbmNlJywKICAgICAgICAgICAgICAgICAgICAgIG15cGFyLmluaXQgPSBOVUxMKQoKCmZpdDAzLmRsbW0kYgpzcXJ0KGZpdDAzLmRsbW0kVikKc3FydChmaXQwMy5kbG1tJHMyKQoKCgp0cnVlX2JldGEgPSBjKGJldGEwLCBiZXRhKQp0cnVlX2JldGEKdHJ1ZV9zaWdtYSA9IGMoc2lnbWFfdSwgc2lnbWFfdl9ob3NwKQp0cnVlX3NpZ21hCmMoc2lnbWFfZSkKCgpwbG90KGZpdDAzLmRsbW0kYiwgdHJ1ZV9iZXRhKQpwb2ludHMoZml0MDMuZGxtbSRiWzFdLCB0cnVlX2JldGFbMV0sIGNvbCA9ICJyZWQiKQphYmxpbmUoYSA9IDAsIGIgPSAxLCBjb2wgPSAiYmx1ZSIsIGx3ZCA9IDIpCgpwbG90KHNxcnQoZml0MDMuZGxtbSRWKSwgdHJ1ZV9zaWdtYSkKcG9pbnRzKHNxcnQoZml0MDMuZGxtbSRWKVsxXSwgdHJ1ZV9zaWdtYVsxXSwgY29sID0gInJlZCIpCmFibGluZShhID0gMCwgYiA9IDEsIGNvbCA9ICJibHVlIiwgbHdkID0gMikKYGBgCgoKCmBgYHtyLCBldmFsID0gRn0KQSA8LSBtYXRyaXgocnVuaWYoMTAyNF4yKSwgbnJvdyA9IDEwMjQpCkIgPC0gbWF0cml4KHJ1bmlmKDEwMjReMiksIG5yb3cgPSAxMDI0KQpzaXplID0gMTAyNAoKYmVuY2g6Om1hcmsobG9nZGV0IDwtIGxvZyhkZXQoZGlhZygxLCBzaXplKSArIHQoQSkgJSolIEEgJSolIEIpKSkKYmVuY2g6Om1hcmsobG9nZGV0IDwtIGFzLm51bWVyaWMoZGV0ZXJtaW5hbnQoZGlhZygxLCBzaXplKSArIHQoQSkgJSolIEEgJSolIEIsIGxvZ2FyaXRobSA9IFRSVUUpJG1vZHVsdXMpKQpgYGAKCgo=